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Abstract 

We combine Differential Evolution, a new technique, with the traditional grid based 
method for optimization of solar neutrino oscillation parameters Am 2 and tan 2 9 for the case 
of two neutrinos. The Differential Evolution is a population based stochastic algorithm for 
optimization of real valued non-linear non-differentiable objective functions that has become 
very popular during the last decade. We calculate well known chi-square (x 2 ) function 
for neutrino oscillations for a grid of the parameters using total event rates of chlorine 
(Homcstake), Gallax+GNO, SAGE, Superkamiokande and SNO detectors and theoretically 
calculated event rates. We find minimum \ 2 values in different regions of the parameter 
space. We explore regions around these minima using Differential Evolution for the fine 
tuning of the parameters allowing even those values of the parameters which do not lie on 
any grid. We note as much as 4 times decrease in x 2 value in the SMA region and even 
better goodness-of-fit as compared to our grid-based results. All this indicates a way out of 
the impasse faced due to CPU limitations of the larger grid method. 



1 Introduction 

The flux of solar neutrino was first measured by Raymond Davis Junior and John N. Bahcall 
at Homestake in late 1960s and a deficit was detected between theory (Standard Solar Model) 
and experiment pQ. This deficit is known as the Solar Neutrino Problem. Several theoretical 
explanations have been given to explain this deficit. One of these is neutrino oscillations, the 
change of electron neutrinos to an other neutrino flavour during their travel from a source point in 
the sun to the detector at the earth surface [2] . There was no experimental proof for the neutrino 
oscillations until 2002 when Sudbury Neutrino Observatory (SNO) provided strong evidence for 
neutrino oscillations [3]. The exact amount of depletion, which may be caused by the neutrino 
oscillations, however, depends upon the neutrino's mass-squared difference Am 2 = m 2 — m 2 (mi 
and m2 being mass eigen-states of two neutrinos) and mixing angle 6, which defines the relation 
between flavor eigen-states and mass eigen-states of the neutrinos, in the interval [0, vr/2]. 

The data from different neutrino experiments have provided the base to explore the field of 
neutrino physics. In the global analysis of solar neutrino data, we calculate theoretically expected 
event rates with oscillations at different detector locations and combine it with experimental 
event rates statistically through the chi-square (x 2 ) function, as defined below by Eq.(TT]) , for a 
grid of values of the parameters Am 2 and tan 2 #. The values of these parameters with minimum 
chi-square in different regions of the parameter space suggest different oscillation solutions. The 
names of these solutions, found in the literature, along with specification of the regions in the 
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parameter space are: Small Mixing Angle (SMA: 10~ 4 < tan 2 9 < 3 x 1(T 2 , 3 x lCT 7 eV 2 < 
Am 2 < 10~ 4 eV 2 ), Large Mixing Angle (LMA: 3 x 1(T 2 < tan 2 6> < 2, 2 x lCT 6 eV 2 < Am 2 < 
10" 3 eV 2 ), Low Probability Low Mass (LOW: 3 x 10~ 2 < tan 2 6> < 2, 10~ 8 eV 2 < Am 2 < 
2 x 10" 6 eV 2 ) and Vacuum Oscillation (VO: 0.1 < tan 2 6» < 1, 10~ n eV 2 < Am 2 < 10" 8 eV 2 ) 
[U . Extensive work has been done on the global analysis of solar neutrino data [SI El El El El 
\TU[ [TT| [T2"l \13\ [T4"l [T5] and now is the era of precision measurement of the neutrino oscillation 
parameters [HI [T7] . 

Traditionally, the whole parameter space (1CT 4 < tan 2 6> < 10, 10~ 13 eV 2 < Am 2 < 10~ 3 eV 2 ) 
is divided into a grid of points by assigning a variable to each parameter and varying its logarithm 
uniformly. The chi-square values are calculated for each point in the parameter space either by 
using 8 B flux constrained by the Standard Solar Model, e.g., BS05(OP) [18] in our case, or by 
using unconstrained 8 B flux [9] where it is varied about the value predicted by the Standard 
Solar Model. The global minimum chi-square value Xmin * s f° un d and 100/3% C.L. (Confidence 
Level) contours are drawn in the tan 2 # — Am 2 plane by joining points with x 2 = Xmin ^X 2 
for different confidence levels. From the chi-square distribution one can easily find that A^ 2 = 
2.28,4.61,5.99,9.21,11.83 for 68%, 90%, 95%, 99% and 99.73% C.L. for two degrees of freedom. 
Minimum chi-square values are found in all the regions and the goodness-of-fit, corresponding 
to each of the minimum chi-square, is calculated. To find the each goodness-of-fit the chi-square 
distribution is used and confidence level 100(1 — /3)%, corresponding to the minimum chi-square 
in the region and the degree of freedom of the analysis, is calculated [U |9]. In our analysis we 
used total event rates of chlorine (Homestake), Gallax+GNO, SAGE, Superkamiokande, SNO 
CC and NC experiments. So the number of degrees of freedom was 4 (6(rates)-2(parameters: 
tan 2 9 and Am 2 )). 

When we use the Differential Evolution (DE), the parameters are randomly selected in the 
given range and checked for a decrease of chisquare, in contrast with the traditional grid based 
method as described in the above paragraph. Thus we selected the vectors with least chi-square 
values, in different regions of the selected grid, as starting points and used DE for the fine tuning 
of the parameters by exploring region around the selected vectors in the parameter space. 

Here in section (2J we define the chi-square (x 2 ) function for the solar neutrino oscillations. 
We use the same x 2 function definition in the algorithm of DE as well as in the traditional 
method. In section [3l we describe algorithm of Differential Evolution along with its salient 
features. In section 2] and [5j we describe results of global analysis by grid and those obtained 
using Differential Evolution respectively. Our conclusions are given in section [6j 

2 Chi-square (x 2 ) Function Definition 

In our x 2 analysis, we used the updated data of total event rates of different solar neutrino 
experiments. We followed the x 2 definition of ref. |19j and included chlorine (Homestake) [20] . 
weighted average of Gallax and GNO [21j, SAGE [22], Superkamiokande [23], SNO CC and SNO 
NC |24j total rates. The expression for the x 2 is given as: 

xLtes= E {Rt-RTWnnV^Rn-RT^ C 1 ) 

J1J2=1,6 

where R 1 - 1 is the theoretically calculated event rate with oscillations at detector j and R^ xp is 
the measured rate. For chlorine, Gallax+GNO and SAGE experiments R th and R ex P are in 
the units of SNU (1 SNU=10~ 36 captures/atom/sec) and for Superkamiokande, SNO CC and 
SNO NC these are used as ratio to SSM Eq.([8]) below. Vj 1 j 2 is the error matrix that contains 
experimental (systematic and statistical) errors and theoretical uncertainties that affect solar 
neutrino fluxes and interaction cross sections. For the calculation of the error matrix Vj 1 j 2 we 
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followed ref. [19] and for updated uncertainties we used ref. [25J. For the calculation of theoretical 
event rates, using Eqs.([3][7]) below, we first found the time average survival probabilities, over the 
whole year, of electron neutrino (P^ e {E u )) (E u is the neutrino energy in MeV) at the detector 
locations for the k th neutrino source and for the grid of 101 x 101 values of ^gr- and tan 2 # 
following the prescriptions described in ref. [9]. For the uniform grid interval distribution we 
used the parameters ^S- and tan 2 # as exponential functions of the variables x\ and X2 as: 



Am 2 
E 

and 



10 (o-i*i-i3) ( 2 ) 



tan 2 = i(r 2 ( 2 - - 025 ^) (3) 

so that discrete values of x\ and X2 from to 100 cover the entire tan 2 # — Am 2 parameter space. 
We used the expression for the average expected event rate in the presence of oscillation in case 
of Chlorine and Gallium detectors given as: 

. dE u \ k {Eu)WAE u ){Pi{E u ))}. (4) 

s=no« J th 

Here E 3 th is the process threshold for the jth detector (j=l,2,3 for Homestake, Gallax+GNO and 
SAGE respectively). The values of energy threshold El, for CI, Ga detectors are 0.814, 0.233 
MeV respectively [26J. <pk ar e the total neutrino fluxes taken from BS05(OP) |18j . For Gallium 
detector all fluxes contribute whereas for Chlorine detector all fluxes except pp flux contribute. 
\k{E u ) are normalized solar neutrino energy spectra for different neutrino sources from the 
sun, taken from refs. \27\ I28j. and a e j is the interaction cross section for v e in the jth detector. 
Numerical data of energy dependent neutrino cross sections for chlorine and gallium experiments 
is available from ref. [27] • Event rates of Chlorine |20j and Gallium |21| 122] experiments and 
those calculated from Eq.@ directly come in the units of SNU. 

Superkamiokande and SNO detectors are sensitive for higher energies, so <^ are the total 8 B 
and hep fluxes for these detectors respectively. The expression of the average expected event 
rate with oscillations for elastic scattering at SK detector is as below: 

Nsk = E <f>* / dE u X k (E u ) x {a ejj (E v )(P* e (E v )) + <J^(E y )[l - (P e k e (E u ))}} . (5) 

k=l,2 J ° 

Here a e and are elastic scattering cross sections for electron and muon neutrinos that we 
took from ref. |29j . 

For the SNO charged-current (CC) reaction, v e d — > e~pp, we calculated event rate using the 
expression: 

Ncc =J2faf dE u \k(E u )a C c(E u ) x (P* e (E u )), (6) 

fe=l,2 ^ 

Here occ is vd CC cross section of which calculational method and updated numerical results 
are given in refs. [30J and |31| respectively. 

The expression for the SNO neutral-current (NC) reaction, v x d — > v x pn {x = e, fi, r), event 
rate is given as: 

K h c = E fa I dE u \ k (E u )a NC {E u ) x ((P* e (E v )) + (P e k a (E u ))). (7) 

fc=l,2 J 

Here a^c is vd NC cross section and {P^ a (E u )} is the time average probability of oscillation 
into any other active neutrino. We used updated version of CC and NC cross section data 
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from the website given in ref. [31]. In case of oscillation of the v e into active neutrino only, 
(P*,(E V )) + (P* a {E u )) = 1 and N$ c is a constant. 

For Superkamiokande [23] and SNO [24] experiments, the event rates come in the unit of 
10 6 em -2 s . We converted these rates into ratios to SSM predicted rate. We also calculated the- 
oretical event rates as ratios to SSM predicted rate in order to cancel out all energy independent 
efficiencies and normalizations [8j. 

th 

= jySSM ( 8 ) 

Here N? SM (j =4,5,6 for SK, SNO CC and SNO NC respectively) is the predicted number 
of events assuming no oscillations. We used the Standard Solar Model BS05(OP) [18] in our 
calculations. Theoretical event rates, so calculated, were used in Eq.([T]) to calculate the chi- 
square function for different points in the tan 2 # — Am 2 parameter space. 



3 Differential Evolution 

Differential Evolution (DE) is a simple population based, stochastic direct search method 
for optimization of real valued, non-linear, non-differentiable objective functions. It was first 
introduced by Storn and Price in 1997 [32] . Differential Evolution proved itself to be the fastest 
evolutionary algorithm when participated in First International IEEE Competition on Evolu- 
tionary Optimization [33J. DE performed better when compared to other optimization methods 
like Annealed Nelder and Mead strategy [34j, Adaptive Simulated Annealing [35], Genetic Al- 
gorithms [36j and Evolution Strategies |37] with regard to number of function evaluations (nfe) 
required to find the global minima. DE algorithm is easy to use, robust and gives consistent 
convergence to the global minimum in consecutive independent trials [32j [38] . 

The general algorithm of DE [39J for minimizing an objective function carries out a number 
of steps. Here we summarize the steps we carried out for minimizing the x 2 function defined 
in section [2J We did optimization of the \ 2 function individually for different regions of the 
parameter space to do one fine tuning in each region. The results of the optimization are 
reported in the section [5] below. 



Step I 

An array of vectors was initialized to define a population of size NP=20 with D=2 parameters as 

x i = x j,i where i = 1, 2, ,NP and j = 1, ..,D. (9) 

The parameters, involved here, are x\ and x-i of Eqs.([2]) and ([3]) on which Am 2 /E and tan 2 # 
depend. Upper and lower bounds {bj,u and individually for different regions of the param- 
eter space described in the introduction section, for the x values were specified and each vector 
i was assigned a value according to 

Xjj = randj(0, 1) • (bjp - hj t £) + b jtL (10) 

where randj € [0, 1] is j th evaluation of a uniform random number generator. The x 2 function 
was calculated for each vector of the population and the vector with least x 2 function value was 
selected as base vector x ro . 



Step II 

Weighted difference of two randomly selected vectors from the population was added to the base 
vector x ro to produce a mutant vector population v, of NP trial vectors. The process is known 
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as mutation. 

Vj = x ro + F ■ (x ri - x r2 ). (11) 

Here the scale factor F £ [0, 2] is a real number that controls the amplification of the differential 
variation. The indices r\,r2 € [1,NP] are randomly chosen integers and are different from r . 

Different variants of DE mutation are denoted by the notation 'DE/x/y/z', where x specifies 
the vector to be mutated which can be "rand" (a randomly chosen vector) or "best" (the vector 
of the lowest \ 2 from the current population), y is the number of difference vectors used and 
z is the crossover scheme. The above mentioned variant Eq. (|lip is DE/best/l/bin, where the 
best member of the current population is perturbed with y=l and the scheme bin indicates 
that the crossover is controlled by a series of independent binomial experiments. The two 
variants, reported in the literature |32|, I38j. very useful for their good convergence properties, 
are DE/rand/l/bin 

Vj = x n + F ■ (x r2 - x r3 ), (12) 

and DE/best/2/bin 

Vj — X ro -)- F • (x ri ~1- X r2 -^T3 ^"7*4 J * 

For our problem, we used the variant DE/best/2/bin Eq. (|13j) for DE mutation, where 2 
difference vectors were added to the base vector. The values of F we used are reported in 
section [5] below. 



Step III 

The parameters of mutant vector population Eq. (jl3p were mixed with the parameters of target 
vectors Eq.([9]) in a process called uniform crossover or discrete recombination. After the cross 
over the trial vector became: 

u . = u .. = l v 3,i If (randj(0, 1) < Cr or j = j rand ), 
* 3 ' % I x j4 otherwise. 

Here Cr € [0, 1] is the cross over probability that controls fraction of the parameters inherited 
from the mutant population (the values of Cr we used are given in section 5), randj G [0, 1] is 
the output of a random number generator and J ran d £ [1)2] is a randomly chosen index. 



Step IV 

The x 2 function was evaluated for each of the trial vector Uj obtained from Eq. (|14p . If the 
trial vector resulted in lower objective function than that of the target vector Xj, it replaced 
the target vector in the following generation. Otherwise the target vector was retained. (This 
operation is called selection.) Thus the target vector for the next generation became: 

A = { Ui If X Y ] - * 2(Xi) ' (15) 

[ Xj otherwise. 

The processes of mutation, crossover and selection were repeated until the optimum was 
achieved or the number of iterations (generations) specified in section [5] were completed. 



4 Analysis from the Selected Grid 

Figure[T]and Table[T]show our best fit oscillation parameters, in different regions, calculated 
using a grid of 101 x 101 points of the parameter space. The symbol of star shows the best 
fit points in the respective regions of the parameter space. Calculations of goodness-of-fit and 
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Solution 


Am 2 




tan 


2 9 


x 2 ■ 

A mm 


g.o.f. 


LMA 


2.512 


• IO" 5 


3.981 • 


io- 1 


0.808 


93.77% 


VAC 


6.31 • 


io- n 


1.00 


10° 


1.268 


86.72% 


LOW 


1.00 


io- 8 


1.122 


•10° 


4.09 


39.46% 


SMA 


6.31 


IO" 6 


1.585 • 


IO" 3 


7.78 


10.01% 



Table 1: Our best-fit values of the oscillation parameters Am 2 , tan 2 8 along with Xmm d.o.f) 
(6(rates)-2(parameters: tan 2 #,Am 2 )) and g.o.f. corresponding to Figure [TJ 




Figure 1: Our global solutions for the total rates. The input data includes total event rates 
of chlorine [20], weighted average of Gallax and GNO (2TJ, SAGE [22], Superkamiokande [23], 
SNO CC and SNO NC [21]. The increasing grey level shows 90%, 95%, 99%, 99.73% C.L. Our 
best-fit points in different regions are marked by stars. 



confidence level are described in the introduction section. We used chi-square function definition 
of section [2j We used 8 B flux constrained by the Standard Solar Model BS05(OP). We saw that 
the point with global minimum or the best fit point in the parameter space lies in the LMA 
region with Am 2 = 2.512 • 10~ 5 eV 2 and tan 2 6 = 3.981 • 10 _1 that is consistent with the results 
found in the literature where SNO data is included in the analysis [101 12] • Before including 
the SNO data the best fit was found in the SMA region [2"6] . 

A selection of a fine grid with larger number of points in the parameter space, of course, will 
give better results but limitations of the CPU time restricted us, like others, to a grid with a 
small number of points. But we point out that even without increasing the number of points in 
the grid we can get lower x 2 an d better g.o.f. by fine tuning of the oscillation parameters using 
DE. We describe what we mean by fine tuning and report our improvements obtained this way 
in the next section. 
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Solution 


Iterations 


Am 2 (eV 2 ) 


tan 2 6» 




Xmin 


xLst 


g.o.f. 




1-7 


2.51189 • 10~ 5 


3.98107 • 


io- 1 


0.808314 




93.77% 




8-9 


2.50999 • 10~ 5 


3.97855 • 


IO" 1 


0.807316 








10-13 


2.40927 • 10~ 5 


3.97684 • 


io- 1 


0.804711 








14 


2.49798 • 10~ 5 


3.97287 • 


io- 1 


0.804564 






LMA 


15 


2.49307 • 10~ 5 


3.97028 • 


io- 1 


0.804289 








16 


2.45884 • 10~ 5 


3.97912 • 


IO" 1 


0.804424 








17-21 


2.46633 • 10~ 5 


3.9751 • 10- 1 


0.803315 








22-29 


2.43264 • 10~ 5 


3.98097 • 


IO" 1 


0.803192 








30-50 


2.45084 • IO" 5 


3.9751 • IO" 1 


0.802953 


0.802953 


93.82% 




1 


6.30957 • HT 11 


1.0 




1.26779 




86.72% 




2-3 


6.64977 • 10~ n 


1 0977 




1.260239 






VAC 


4-6 


6.70641 • 10~ n 


1.01912 




1.25977 








7-43 


6.6423 • 10~ n 


0.993134 




1.25948 








44-45 


6.6423 • 10~ n 


0.998353 




1.25945 








46-50 


6.70041 • 10~ n 


0.99326 




1.25939 


1.25939 


86.82% 




1-3 


6.30957 • 10~ e 


1.58489 • 


IO -3 


7.77933 




10.01% 


SMA 


4 


6.19532 • 10~ 6 


1.64599 • 


10~ 3 


6.24974 








5 


6.10563 • 10~ 6 


1.48276 • 


10~ 3 


5.89341 








6-50 


5.48095 • 10~ 6 


1.72371 • 


10~ 3 


1.86456 


1.86456 


75.97% 




1 


1.0 • 10~ 8 


1.12208 




4.18897 




39.46% 




2-4 


2.37807 • 10~ 8 


1.03198 




3.98339 








5-9 


2.95404 • 10~ 8 


1.03198 




3.97624 






LOW 


10 


3.3042 • 10~ 8 


1.03069 




3.97605 








11 


2.80796 • 10~ 8 


1.02741 




3.9728 








12-20 


3.17357 • 10~ 8 


1.02741 




3.96267 








21-50 


3.14543 • 10~ 8 


1.02723 




3.96125 


3.96125 


41.12% 



Table 2: The results of the oscillation parameters during different iterations of the DE algorithm. 
The improved values of the oscillation parameters Am 2 , tan 2 9 along with Xbest d.o.f) and g.o.f. 
using Differential Evolution strategy DE/best/2/bin corresponding to Figure [2] are presented. 
Note in the 1 st row of different regions, the points with minimum chi-square given in table Q] are 
taken as first members of the population arrays. The other members of the arrays, for different 
regions, are selected randomly using DE. 

5 Optimization of the Chi-square Function using DE 

We have described algorithm of the Differential Evolution in detail in section [3l We wrote 
the subroutine of the chi-square function, denoted by x 2 > following the definition of chi-square 
in section [2j that depends on xi and %2 and used it as objective function of the DE algorithm. 
We combined the traditional grid-based method with DE in two aspects: First, we used the 
survival probabilities {P^ e (E u )) already calculated for the discrete values of x\ and X2 for our 
grid of 101 x 101 points of the parameter space and interpolated them to the continuous values 
of x\ and X2 to calculate event rates and chi-square function in DE algorithm. We used cubic 
polynomial fit for the interpolation purpose to fit the data. Second, we used the points with 
minimum chi-square in different regions of the selected grid Table [1] as the starting points (and 
members of the respective population array) and explored the space around them for the fine 
tuning. That is, we searched for the points with smaller x 2 values and better goodness-of-fit of 
the oscillation parameters. 
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0.0000255 




0.0000235 L 

0.36 0.37 0.38 0.39 0.40 0.41 0.42 0.43 



Taa*0 



7.x 10 




5.xlO" sL 

0.0010 0.0015 0.0020 0.0025 0.0030 

Tan 1 e 




0.95 1.00 1.05 1.10 1.15 
Tan J 



Ox 10 




1.10 



Tan' 9 



Figure 2: Track of the DE algorithm for optima in different regions using the strategy 
DE/best/2/bin. The square symbol shows the best point of the 101 x 101 grid and trian- 
gle symbol shows the best point after fine tuning using DE. 



In our analysis, the values of DE control variables F and CR were taken as 0.4 and 0.9 
respectively for the LMA, SMA and VAC regions. For the LOW region F and CR were both 
taken as 0.3 for better convergence. Maximum number of iterations were taken to be 50 for all 
regions. We took the best point in a region of the 101 x 101 grid in Table [T] as the first member 
of the population in the first iteration and used the strategy DE/best/2/bin for DE mutation 
in all the remaining iterations/generations. The steps of DE algorithm, described in section [3l 
are repeated for the number of iterations specified. 

Table [2] and Figure [2] show the results in different regions during and after fine tuning of the 
oscillation parameters using Differential Evolution. The value of Xmm persisted, rejecting all the 
mutations, for the iterations mentioned in column 2 of Table [2j Accepted mutations resulted 
in new vectors whose components are given in column 3 and 4 of the following rows. Xbest * s 
the minimum chi-square value we obtained in the region specified. In comparison to the results 
of Table [IJ we note here as much as 4 times decrease in the Xmin °f the SMA region after fine 
tuning using DE along with a small decrease in all the other regions. Different vectors in Figure 
[2] show the track of DE algorithm for optima in different regions during iterations specified in 
Table H 
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6 Conclusions 



Fine tuning of the neutrino oscillation parameters using Differential Evolution has been 
introduced as a solution to the impasse faced due to CPU limitations of the larger grid alterna- 
tive. We can explore the parameter space deeply due to real nature of the parameters X\ and X2 
using DE in contrast to discrete nature of these parameters in the traditional grid based method. 
We conclude that combination of Differential Evolution along with traditional method provides 
smaller chi-square values and better goodness-of-fit of the neutrino oscillation parameters in 
different regions of the parameter space. We also note a significant change in the results of Xmin 
and g.o.f. in the SMA region after the fine tuning using DE. Even though it is a local decrease, 
it indicates importance of the exploration of the points within the grid and the efficiency that 
can be achieved through DE. 
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